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An efficient and accurate quantum algorithm for the Dirac equation 

Jeffrey Yepez^ 
Air Force Research Laboratory 
29 Randolph Road, Hanscom Field, Massachusetts 01731 
(Dated: March 25, 2002) 

An efficient quantum algorithm for the many-body three-dimensional Dirac equation is presented. 
Its computational complexity is dominantly linear in the number of qubits used to spatially resolve 
the 4-spinor wave function. 
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The first quantum algorithm to compute a path inte- 
gral was introduced by Feynman in 1965. At that time 
he did not refer to it as such since it was not until in 1982 
that he proposed the idea of using a quantum computer 
to efficiently simulate quantum mechanical many-body 
dynamics ||l|. In the second chapter of his manuscript 
on path integrals published with Hibbs the problem 
is given to prove that the one-dimensional (ID) Dirac 
equation can be modeled by summing over all the pos- 
sible zigzag paths of a single-speed particle traveling at 
the speed of light in a discrete two-dimensional space- 
time hopping from lattice site to lattice site and flipping 
left or right according to a deterministic local interaction 
rule. The amplitude a particular path contributes to the 
kernel is proportional to the number of its "collisions" or 
corners. In this way, the time evolution of the 2-spinor 
field of a single quantum particle is modeled by a "gas" of 
particles computing all paths simultaneously in a time- 
explicit fashion. This discrete and parallelized process 
of local collisions and the streaming along the lattice is 
described by a continuous effective field theory, the ID 
Dirac equation, in the limit of the infinite lattice reso- 
lution. A solution to Feynman's "quantum lattice gas" 
problem was published in 1984 by Jacobson and Schul- 
man (|. Here we give a three-dimensional (3D) solution. 

In 1994, Bialynicki-Birula proposed a discrete model of 
the 3D Dirac equation implemented on a body-centered 
cubic lattice Q. However, this model is Ist-order con- 
vergent (doubling the grid resolution merely doubles the 
numerical accuracy) , problematic when modeling particle 
dynamics in an external potential. Although the model 
is unitary, it is specified using non-unitary matrices and 
requires ad hoc lattice partitioning if implemented in par- 
allel. Furthermore, Bialynicki-Birula addresses only the 
one-body problem. Meyer published a subequent series of 
papers on the ID quantum lattice gas algorithm, equiv- 
alent to Bialynicki-Birula's algorithm, and cast in the 
form of Feynman's original model jsj . Meyer contributed 
ID one-body numerical simulations and addressed the 
non-interacting lattice or checker-board problem using 
an additional rest particle. Yet he too did not address 



the many-body case nor the low-order numerical conver- 
gence issue. 

Contemporaneously with Meyer, Succi published a se- 
ries of papers on this subject emphasing the analogy be- 
tween quantum mechanics and fluid mechanics: the con- 
nection between the Dirac equation and the Schroedinger 
equation to that between the kinetic Boltzmann equa- 
tion and the Navier-Stokes equation of hydrodynamics 
Succi's quantum lattice gas model on a cubic lattice 
for the 3D Dirac equation has, at the "kinetic" level, the 
particles undergoing mixing during free propagation and 
is again similar to Bialynicki-Birula's model. Succi dis- 
cusses the many-body case, but his algorithm runs into 
an "exponential complexity wall" . 

Our quantum lattice-gas algorithm for the 3D Dirac 
equation is suited to direct implementation on a quantum 
computer using only two-qubit quantum gates and effi- 
ciently handles the many-body problem. For pedagoical 
purposes, first we state the simplest quantum lattice- 
gas algorithm on a cubic lattice. Then, we introduce 
an improved version that remedies two difficiencies: the 
checkboard problem of non-interacting sublattices and 
the low-order convergence. Finally, we recast our quan- 
tum algorithm to handle the many-body case in a second- 
quantized representation. 

The relativistic quantum mechanical wave equation for 
a free particle is the linear Dirac equation 



c > aidiip -I- 1— —pip, 



(1) 



where ip is a 4-spinor and the matrices ai and P sat- 
isfy the constraints af — 1, P'^ — 1, {ai,aj} ~ 0, and 
{/3, ai\ = so that (|l]) is equivalent to the Klein-Gordon 
wave equation. Since the 2x2 Pauli matrices Ui for 
i = x,y,z satisfy these constraints, we can express at 
and 13 as tensor products ai = a ® Ui and (3 = b ® 1, 
where a, b, can be any two different Pauli matrices and 
where 1 is the 2x2 identity matrix. We choose a = 
and b = ax- Then, the Dirac equation is 
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(2) 
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With the wave function defined on an infinite resolution 
cubical lattice at times separated by an infinitesimal du- 
ration St with the grid cell size the infinitesimal length 
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Sr = c5t, the Heisenberg evolution 

V^' = V + (5V = e^''"^®'""'''^'"*'^^'^*'"^®^ (3) 

corresponds exactly to (0) in the relativistic limit where 
fvjj mc^ and hk ~ mc. 
The 4x4 matrix 
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(4) 



operating with the z spatial derivative in (|^) is diagonal 
whereas the matrices az ® Ux and CTz ® Uy for the x and 
y partial derivatives, respectively, are not diagonal. 

We would like to transform (^) in such a way that all 
the matrices operating with the spatial partial deriva- 
tives are diagonal (and hence correspond to infinitesimal 
shifting along the orthogonal lattice directions). To do 
this, we need the two identities: 



(5) 



that follow from e^J'^' = (1 + ^fj) provided e is i 



finitesimal. Then, using the identity 1 ® e = e" 
the 2-spinor similarity transformations can be gener- 
alized to 4-spinor transformatons 



(1- 



(6) 
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which we will use to diagonalize the x and y spatial 
derivative operators in (0). Using m) and defining 
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(7) 



(8) 



(9) 



and Si 



the spatial displacement operators 
in the Heisenberg representation of the evolution equa- 
tion (|^) can be written 
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w y 2L 



(10) 



so the evolution equation itself can be rewritten as 

V,' = Yi'^SxYP'x^^^'SyX^^^S.X^^U. (11) 



This has the form of a quantum lattice-gas algorithm 

(2) (2) 

with local interaction ( "collision" ) operators , Y^ ' 



and X^^ll^^ , as well as lattice-directed displacement 

K 

("streaming") operators Sx, Sy, and 5*2. 

For numerical purposes, we would like to represent the 
wave function on a finite resolution grid with cell size 
Sr — > Ar and update time 6t — s- At. In this approxi- 
mation, the wave function becomes a discrete field that 
exists only at the spacetime grid points xe for £ = 1 . . . L'^ 
and tn for n = 0, 1, 2, . . . 



i>ixe,tn) 



/a{xi,t„ 

f3{xi,tn) 

fj.{xe,tn) 

\v{xi>,tn)j 



(12) 



and the operators Si for i — x^y, ov z induce a finite 
displacement Si'4'{x) 'ip{xe + ® CTz Arii) of the com- 
ponents of the 4-spinor only along lattice directions: 



Sx'4'ixe,yi,ze) 
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(13) 



and 



Sz'>P{xe,yi,Zi) 



fa{xe,yi + Ar, zg) 
I3{xi,yi - Ar,Zi) 
H{xi,yi - Ar,Zi) 

\v{xi>,yi + Ar,Zi) 



fa{xi,ye,ze + Ar)\ 
P(xi,yi,Zi - Ar) 
^i{xi,yi,Zi - Ar) 

\i^{xe,ye.,ze + Ar)/ 



(14) 



(15) 



These streaming operators are classical data shifting op- 
erators causing global permutations of the components 
of the 4-spinor wave function across the lattice and on 
a quantum computer can be implemented by 2-qubit lo- 
cal swap operators The collision operators act in- 
dependently on each node of the lattice and cause local 
quantum entanglement between component pairs of the 
4-spinor. The streaming operators in turn propagate this 
local on-site entanglement to next nearest neighbors so 
that eventually quantum entanglement covers the entire 
lattice. 

It is possible to rewrite ( [III ) as a finite difference equa- 
tion on a body-centered cubical lattice. The resulting 
set of coupled finite difference equations are similar to 
the finite difference representation of the 3D Dirac equa- 
tion given by Bialynicki-Birula in 1994 A drawback of 
expressing the algorithm as a finite-difference equation is 
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its unsuitability for a quantum computer implementation 
using two-qubit quantum gates whereas our manifestly 
unitary expression ( |ll| ) is suitable. 

A continuous effective field theory for ip — (a, (3, fi, u) 
follows in the continuum limit of the emergent finite- 
difference equations by Taylor expanding in Ar = x^+i — 
xi and in = — t„. We obtain 



and 



c 




At 



(16) 



— r, +0(cAr,Ai), 



which is exactly the Dirac equation (|l|) when Ai ^ Ar ^ 
e are infinitesimal and when the partial derivative with 
respect to time is defined as dtij} = ^ ■ (11) gives 



rise to perfectly unitary evolution of the discretized wave 
function and, therefore, is an unconditionally stable nu- 
merical algorithm. The effective field theory (|l^) is f st- 
order convergent in space. 

Our basic approach to improve the accuracy of the 
quantum algorithm is to set the grid size Ar to be smaller 
than the Compton wavelength A = ^ of the modeled 
particle 



Ar 



(17) 



and to introduce a small temporal scale that is much 
smaller than - 



Ai 



(18) 



The diffusive ordering condition of spatial and tempo- 
ral fluctuations typical of random walk processes, Ar^ = 
vAt^ provides a context to understand the scaling be- 
havior of the small parameter e. According to (17) and 
( p^ ) , the diffusive transport coefficient is v = and the 
particle velocity is = | , which approaches infinity as 
e — > 0. In this limit, the velocity of the modeled quan- 
tum particle is relatively small, hence the resulting effec- 
tive field theory should correspond to the non-relativistic 
limit of the Dirac equation as e ^ 0. 

To diagonalize the streaming operators in ( p^ , we used 
a fixed and finite rotation angle ^ independent of the 
grid resolution. We will now diagonalize the streaming 
operators using a small rotation angle proportional to 
At. By (|l|), the rotation angle is 6* = ^^^^ = which 
is dependent on the grid resolution. The displacement 
operators in the Dirac equation (|l]) can be represented 
by interleaving streaming and collision operators on a 
cubical lattice as follows: 



F = 02,4^^(2) 02.4^^(2)1 c.1,3ta(2) ol,3^(2)t 

J-^x — X — X — X — X — 



(19) 



z^aySrdy 



x-(2) Ql,3y(2)tQl,3y(2) 

2^2 y 2 



(20) 

where the superscripts on the streaming operators refer 
to individual components of the 4-spinor. The streaming 
operators Si = S^'^S^'* in (13) and (|lj) are now sep- 
arated by collision operators. This kind of interleaving 
of streaming and collision operators removes the spuri- 
ous check-board effect of independent sublattice dynam- 
ics that otherwise occurs fl P ho . So far we treated the 



but 



non-diagonal operators e"^'^"^^ and e'^"^'^^^^^^ 
not the displacement operator ^'^^^'^zSrd^ because no such 
improvement exists since it is diagonal. However, if in- 
stead of using the Dirac matrix az ® <Jzi we use an al- 
ternative non-diagonal representation for the z-direction 
partial derivative, then we can again employ interleaving. 
Therefore, we consider this alternate form of the Dirac 
equation 

dt^i = caz(^(Jxdx'ip+caz(S>crydyip+cay'^ldz4'+i<^x'9^^—''P- 

. 

Now the displacement operator in ( |2l| ) for the z-direction 
can be re-expressed in a fashion similar to (|l^) and ( [20| ) 

as 

^cTyf^urd, _^ = s'^'^x9''' s^'^x^^^^ s^'^x''^^ s^'''^x9'^^ 

(22) 

Then instead of (^, we use (|l|), and (H) for an 

improved quantum algorithm 



^{t + At)^ExEyEzXP^^P{t). 



(23) 



to 



In ( p^ ) we have appended a collision operator Xe^^^ 
produce the mass term in the Dirac equation. 

It is possible to derive a finite-difference equation rep- 
resentation of the quantum lattice-gas algorithm ( p3|) by 
carrying out all the collision and streaming operations 
symbolically. The result is no longer expressible on the 
body-centered cubic lattice. Nevertheless, once again, a 
continuous effective field theory for ijj = [a, (3, ^,v) fol- 
lows in the continuum limit and Taylor expanding in Ar 
and in At: 



Ard. 



+ lArdz 
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—a 



(24) 



0(Ar^A^), 



which approximates the Dirac equation (^TJ) when At is 
small. The effective field theory ( p^ is Ist-order conver- 
gent because of the error term 0{At). With the evolu- 
tion operator E = E^EyEz, we define the dual operator 
E = El^^E^_yE^__^, by taking the adjoint of the collision 
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FIG. 1: L2 norm error ^ iT,f=iWi^i)\'^ 
versus grid resolution 5x ■ 

sizes from L = 8 to L = 32768. The error curve's slope of the 
original and improved algorithm is 0.5 (dashed line) and 2.5 (solid 
line), respectively. This demonstrates the high numerical accuracy 
of the improved quantum algorithm. 



operators and reversing the streaming directions. Then, 
it is possible use a symmetrized evolution operator [nOl 



ip{t + At) ^ EEe-^^ tPit) 



(25) 



which is better than 2nd-order accurate in space, as 
demonstrated in Fig. |[ 

Our quantum algorithm for the many-body Dirac 
equation can be expressed in terms of 2-qubit gates 
that conserve particle number acting on an ini- 
tial ket with 4 qubits per lattice node, j^I') — 

^r=i \Qiir))\q2{r))\q3{r))\qi{r)). With 4, a<„ and h = 
a^Ua denoting the creation, annihilation, and number 
operator, respectively, of the ath qubit (1 < a < 4i^), 



the collision operators are 

Xafi — 1 ~ ism 9{al^ap + d pQa) + {cos 9 ~ l){na + fifs) 



'2 cos Oflaflp 



(26) 



Yap = 1 + s\n6{a)^dj3 — a\aa) + {cosd — l){na + fip) 



-2 cos Ohaflll, 



(27) 



where a and (3 index different qubits at the same site. 
Then are rewritten as X^^^ ^ ^13^24, ' ^ 

X12X34, and yj^-* Y12Y34,. Hence, 2L'^ applications of 
either Xap or Yap are required for a single collision step. 
Streaming occurs by successive application of the inter- 
change operator 5^1^ = 1-1- aJjOiy + aj,a^ — + fii^ [||. 

{L — 1)^ number of applications of 5^1/ (/i refers to one 
qubit-component at some site and v to the same compo- 
nent at its neighboring site) are required to stream one 
qubit-component along a cubic lattice direction. The to- 
tal evolution operator E is the product of collision oper- 
ators X and Y and streaming operators S corresponding 
to algorithm (11) or some variant of (|2^) depending on 
the desired degree of numerical accuracy. With the new 
ket \'i''{t + At)) ^ E\^{t)), the resulting probability of 
finding a particle at site x is P(x) — 
where a the index of the 1st qubit at x. 

The computational complexity of one time step scales 
as C = Pc^L^ + ps{L — 1)^, where pc and ps are the num- 
ber of operations per node for collisions and streaming. 
For the simplest algorithm ( [ll|), pc — 5 and ps — 12, and 
for the improved algorithm (ESh, pc = 13 and ps = 24. 
Both Pc and ps double when we use a symmetrized rule 
like (^) but are a fixed-cost overhead. With Q = AL^ 
qubits, the size of the Hilbert space is exponential 2^5, 
whereas the complexity C = + Ps[Q — |(2(5)3 -|- 
|(2(3)^ — 1] for all the versions of our quantum algo- 
rithm is dominantly linear in Q. 

I would like to thank the Air Force Office of Scientific 
Research for supporting this work. 
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